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We present dynamical transport calculations based on a tight-binding approximation to adiabatic time-dependent 
density functional theory (TD-DFTB). The reduced device density matrix is propagated through the Liouville-von 
Neumann equation. For the model system, 1,4-benzenediol coupled to aluminum leads, we are able to confirm the 
equality of the steady state current resulting from a time-dependent calculation to a static calculation in the conven- 
tional Landauer framework. We also investigate the response of the junction subjected to alternating bias voltages 
with frequencies up to the optical regime. Here we can clearly identify capacitive behaviour of the molecular device 
and a significant resonant enhancement of the conductance. The results are interpreted using an analytical single level 
model comparing the device transmission and admittance. In order to aid future calculations under alternating bias, 
we shortly review the use of Fourier transform techniques to obtain the full frequency response of the device from a 
single current trace. 



1 Introduction The field of quantum transport at the 
molecular scale significantly diversified over the last years 
|[Tl|2l|3]. While the interest was initially to measure the con- 
ductance across individual molecules in an accurate and 
reproducible fashion, current topics involve spin transport 
ID, molecular transistors Q, thermoelectric effects |6,7| 
and device heating |[8]|9l. On the theoretical side much 
progress was achieved using Green's function methods in 
the energy domain |10|. Time domain methods, on the 
contrary, promise easy access to dynamical properties, like 
ac transport, light-induced effects and higher harmonics in 
the current fTT|. In this contribution we report on results 
of such a method based on approximate time-dependent 
density functional theory, termed TD-DFTB |12 13|. The 
scheme allows to perform dynamical transport simulations 
of realistic devices taking the electronic structure of molecule 
and leads into full account. Extending an earlier study on 
a similar topic [ 14 1, we first ask the question whether time 
and energy domain methods provide the same answer for 
the steady state dc current. We continue with a discussion 
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of alternate currents and focus here especially on resonant 
enhancement of the admittance beyond the low frequency 
regime commonly studied. 

2 Method In the following we present a brief descrip- 
tion of our simulation method. A more detailed derivation 
and justification of the present scheme may be found in the 
original articles |15| and |13|. We assume a setup of the 
molecular electronic device as depicted in Fig.[T] The peri- 
odic left (L) and right (R) lead extend to infinity and are in 
thermal equilibrium at the chemical potential iJia=L,R with 
/iL = M-R at t = 0. At i > 0, a time-dependent bias po- 
tential V{t) is applied that drives the central device region 
(D) out of equilibrium and leads to a time-dependent cur- 
rent. Instead of working with the full infinite system, one 
can derive a Liouville equation for the device region only 
IITSl (in atomic units): 






[Hit),a{t)] 



a=L,R 



(1) 
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Here a{t) denotes the one-particle density matrix for 
the device region in a basis of localized atom-centered ba- 
sis functions 4>fj,{r). The Hamiltonian H{t) is given in the 
adiabatic approximation of TDDFT 1161 and depends on 
the electron density p{r,t), while Qa{t) incorporates all 
effects due to the metallic leads, especially also dephasing 
and dissipation. Numerically tractable and explicit forms 
for this term can be obtained from non-equilibriums Green's 
function theory in the wide band limit (WBL). As shown 
by Zheng et al. 1 15], Qa then takes the form: 

g„(t) = iK,CT(t)] + {r„,a(i)} + ir„(i), (2) 

where Aa describes the change of the device energy lev- 
els due to the presence of lead a, while r^-t renders the 
lifetime of the molecular levels finiteO] Both matrices are 
evaluated from first principles and depend on the device- 
lead interaction and the lead surface density of states. The 
term Ka involves only known quantities besides the time- 
dependent bias potential V{t) and hence Eq.fTI represents 
a closed equation that can be numerically integrated by 
conventional Runge-Kutta methods. To this end, the initial 
density matrix at t = may be obtained without further 
approximations from equilibrium Green's function the- 
ory in the WBL. As also shown in reference 1 15|, knowl- 
edge of Qa {t) allows one to compute the time-dependent 
particle current /„ (i) through the left or right device-lead 
interface according to 



/„(t) = -Tr[Q,(i)]. 



(3) 



In practical simulations the time step has to be chosen 
in the attosecond regime in order to resolve the electron 
dynamics accurately. This limits the accessible device di- 
mensions and total simulation time significantly. We there- 
fore adapted the scheme from above to the time-dependent 
density functional based tight-binding (TD-DFTB) method 
fll7lll2l[T3l . In essence, the TDDFT Hamiltonian matrix 
H^v{t) is replaced by 

H^u{t) = {<P^\H[po]\<P,) (4) 

The first term on the right hand side is the DFT Hamil- 
tonian evaluated at a time independent reference density 
po = X^A PA, taken to be a sum of atomic densities pA for 
each atom in the device region. These densities and hence 
also the matrix elements can be computed beforehand. The 
second term involves the overlap 5*^^ of the basis functions 
and takes the deviation of the electrostatic potential from 
the reference into account. The potential Va on the atoms 
in the device region is computed at each time step from a 
Poisson equation with boundary conditions determined by 
the given bias potential in the leads. The charge density 



1^ 



Molecular Device 

Sl Sr 










' With respect to the article by Zheng et al. fTsl, the designa- 
tion of Fa and Aa is interchanged here. Square and curly brack- 
ets indicate a commutator and anti-commutator, respectively. 



Figure 1 Schematic setup of the molecular device shown 
together with our test system. Only the atoms in the device 
region are shown. 



required in this process is computed from the density ma- 
trix a{t) |13|. Besides this adaption in the Hamiltonian, 
we follow the formalism of Zheng et al. without further 
modifications. 

3 Results 

3.1 Approach to steady state We applied the TD- 
DFTB scheme to the junction depicted in Fig.[T] The 1,4- 
benzenediol molecule was optimized with passivating hy- 
drogens in vacuum at the DFTB level and then symmet- 
rically positioned inbetween Al nanowires of finite cross 
sections. The device region consists of the molecule and 36 
additional Al atoms, while the simulation cell for the leads 
included 72 atoms. The latter is periodically replicated to 
+00 and — oo for the right and left lead, respectively, in or- 
der to compute the surface Green's function and WBL pa- 
rameters (see Eq. l2| at zero bias. The basis set is given by 
one s-type atomic orbital for H and one s-type and three 
p-type orbitals for the other elements. The Perdew-Burke- 
Ernzerhof exchange-correlation functional |18| is used in 
all calculations. This model structure was already used in 
(TJI as well as in the first principles TDDFT study 1 14|, so 
that benchmark data is available for comparison. Transport 
through benzenediol is typical for conjugated molecules in 
many respects. The transmission at the Fermi energy Ep 
is rather small (T « 0.01) and transport occurs through the 
tails of the tt and tt* frontier orbitals. 

In Fig. |2] we plot the time-dependent current through 
the left and right molecule-lead interface. Here and in the 
following we integrate Eq.fTlwith a time step of 2 as using 
a 4-th order Runge-Kutta method. The bias voltage of 3.5 
V is applied to the left lead only and turned on exponen- 
tially with a time constant of 0.5 fs. The current initially 
overshoots, oscillates and settles into the steady state only 
after several fs, long after the bias potential nearly reached 
its maximum. Earlier we have shown |13|, that the initial 
transients depend on the time constant of the exponential 
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Figure 2 Absolute value of the time-dependent current 
through the left (/l) and right (Ij^) interface of the molec- 
ular junction depicted in Fig. [T] The inset shows the bias 
potential V{t) = Vb[l - cxp(-i/T)] with Vq = 3.5 V and 
T = 0.5 fs. 



turn on, but not the asymptotic value of the current. In 
addition, we could relate the decay time of the oscillations 
to the imaginary part of the self energy of the device. Well 
coupled junctions reach the steady state earlier, whereas 
weakly coupled junctions feature persistent oscillations 
(see also |19|). As can also be seen in Fig.l2j the absolute 
values of the currents through the left and right interface 
equal each other asymptotically, but not in the transient 
phase of the simulation. Indeed, the particle current is a 
conserved quantity only in the dc limit. Under ac driving 
the device may become charged and one has to consider 
both particle current and displacement current f20l. By 
monitoring the total device charge as a function of time, 
we verified that the latter indeed compensates for the dif- 
ference between |/L(i)| and |/ii;(t)|. 

An interesting question is now, whether the asymptotic 
current /^ — limt_j.oo I{t) from the time-dependent sim- 
ulation equals the current obtained from a conventional 
static calculation in the Landauer formalism. In the latter 
approach the current is given by the energy integral 



/ = Go / dE [f{E, fiL) - f{E, Mfl)] T{E, V) 

T{E, V) = Tr [G^rnG'^rL] , (5) 

with f{E,^) denoting Fermi distribution functions with 

IJ-L—IJ-R = V^, the quantum of conductance Go ~ 77.48/LtS, 
and the bias dependent transmission function T{E, V) fTOl. 
The retarded (G^) and advanced (G") device Green's func- 
tions depend on the Hamiltonian and charge density n(r). 
Since n(r) depends itself on G^ as well as on the applied 
bias, a self-consistent determination of all quantities is re- 
quired. It is not a-priori evident, that the currents given by 



Figure 3 Asymptotic time-dependent current ij^) and 
current in the Landauer formalism (/negf) for 1,4- 
benzenediol as a function of applied bias. The values for 
/^ have been obtained from simulations with a total sim- 
ulation time tmax of 20 fs and a bias potential V{t) ~ 
Vo[l — cxp(— t/r)] with T= 0.5 fs. The current has been 
averaged over the last 2 fs. The wide band approximation 
was also employed in the Landauer calculations. The line 
is a guide to the eye. 



Eq. [3] and Eq. l5] are identical. We have recently discussed 
this question in great detail in the context of first-principles 
TDDFT |14|. Here we perform similar simulations using 
the TD-DFTB method in order to show that our findings 
are not restricted to a specific choice of the Hamiltonian. 
In Fig. [3] we compare the asymptotic currents /^ from 
several time-dependent simulations at different bias values 
with the corresponding values from Eq. [5] Despite signif- 
icant formal and also algorithmic differences between the 
two approaches, one can observe nearly identical values 
over the full bias range. Like in Ref. 1 14|, we conclude that 
time-dependent simulations do in general not offer addi- 
tional or more accurate information when the interest is in 
steady state properties^ As we discuss in the next section, 
there is however an important computational advantage for 
ac transport. 

3.2 Admittance from current traces Starting with 
the work of Fu and Dudley ll24l . several studies addressed 
the response of meso- and nanoscopic devices to an alter- 
nating bias potential 1 25, 26, 27 1 . In recent years approaches 
based on energy domain Green's functions became espe- 
cially popular II28II29|[301[3T1 . but also time domain tech- 
niques, as presented here, allow for the efficient evaluation 
of the admittance ll32l. 



" This statement holds for conventional local and semi-local 
functionals of the density. For non-local functionals differences 
with respect to the Landauer approach have been predicted 1211 
I22I23I . 
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To this end, the Fourier transforrrj^of bias and current 
is numerically evaluated, e.g.. 



V{u:) 



V{t) exp{iujt)dt 



(6) 



to yield the complex admittance Y{uj) = I{uj)/V{uj). In 
electronic circuit theory, the real and imaginary parts of Y 
are also often termed conductance (G ~ Re{Y)) and su- 
ceptance {B = Im(y)), respectively. With the choice for 
the sign of the Fourier transform from above (Eq. l6]l, ca- 
pacitive devices feature a negative susceptance, while in- 
ductive behaviour is characterized by positive values of B. 
We applied this approach to the 1,4-benzenediol junc- 
tion and experimented with different choices for the tem- 
poral profile of the bias potential. In principle, the form of 
V{t) is arbitrary as long as the amplitude is small enough 
to remain in the linear response regime and the support of 
its Fourier transform is sufficiently large. Fig. |4]depicts the 
absolute value of the admittance |i^(w)| for different func- 
tions V{t). As reference, we perform simulations with a 
harmonic bias V{t) = Vq sin(a;t) for different discrete val- 
ues of u and determine the amplitude of I{t) after the ini- 
tial transients have died out. A sample simulation is shown 
in Fig. Is] Inspection of Fig. [4] re veals that an exponential 
turn-on of the form V{t) — Vb[l — exp(— t/T)] provides 
a reasonable estimate for the general features in the ad- 
mittance, but fails to convince on a quantitative level. The 
reason is that the Fourier transform does not exist in the 
limit u! —>■ 0, unless one artificially damps V{t) by a factor 
exp{—rt) to enforce convergence. For small values of F 
the admittance differs strongly from the reference, while 
for larger values the dc limit is overestimated. In response 
calculations for optical properties one often uses a Dirac 
delta function or Lorentzian as a perturbation. Here, the 
Fourier transform exists for all ut and no artificial broaden- 
ing is required. Results for the bias potential 



V{t) 



Vn 



7 



TT (i-to)2+72' 



(7) 



show excellent agreement with the reference data over nearly 
the full frequency range, especially also in the dc limit. A 
benefit with respect to the simulations at discrete frequen- 
cies is that the Fourier transform technique requires only a 
single run to evaluate the full admittance. In the following 
we therefore continue with this choice. 

After this more technical discussion we now analyze 
the admittance in more detail. Fig. l6] a) shows the con- 
ductance and susceptance of 1,4-benzenediol. For small 
frequencies, the negative values of the latter indicate ca- 
pacitive behaviour of the junction. This is in line with the 
simulations shown in Fig. l5] where the current leads the 
voltage signal. The negative susceptance can be rational- 
ized by inspection of the transmission T(E,0) (Fig.l6]c)) of 

' Since V{t) — for f < 0, this is equivalent to the Laplace 
transform with imaginary argument. 




Reference 

Lorentzian 

Exp. with r = 0.1 fs"' 

Exp. with r = 0.01 fs"' 
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Figure 4 Absolute value of admittance |F(w)| in units of 
Go as a function of frequency in units of [eW/h]. Results 
are given for the harmonic perturbation with discrete fre- 
quencies (Reference) and using Fourier transforms with 
exponential form and damping of the bias (Vq = 1 mV, 
T = 0.2 fs, imax = 50 fs) as well as with Lorentzian form 
(Vb = 0.1 mV, 7 = 0.2 fs, to ^2 fs, i^ax = 50 fs). 




20 30 
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Figure 5 time-dependent current due to the ac bias V{t) 
Vo sin{Lot) (Vo = 0.5 fiV, oj = 0.59 eW/h, i„ax = 50 fs). 



the junction. For small frequencies, only the region around 
the Fermi energy is relevant in the linear response regime. 
Here the transmission is low and the current effectively 
blocked, similar to a macroscopic capacitor In a classical 
RC circuit, the admittance is given by 



Y'^^iu) = -iujC + u^C^R 



(8) 



up to second order in the frequency 1281 . As seen in Fig.|6] 
a), real and imaginary part of Y{u!) show for small fre- 
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Frequency [eV {2n/h)] 



Figure 6 a) Real and imaginary part of the admittance 
Y{uj) for 1,4-benzenediol. b) Analytical results for the 
one-level model of Fu and Dudley |24| with parameters 
AE — 2 eW and 7 = 0.25 eV. c) Transmission T(E,0) of 
1 ,4-benzenediol. 



for a single channel. Still a rough assignment of the res- 
onances to individual molecular states becomes possible. 
The first resonance likely originates from the unoccupied 
states 1.5 eV above and the occupied states 2.2 eV below 
the Fermi energy. According to the analytical results, states 
with smaller broadening contribute less to the admittance. 
The resonance around 4 eV is therefore assigned to the 
broad resonance in the transmission around this energy. 
The underlying physical picture is that the time-dependent 
bias potential opens new transport channels at E+Tilu and 
E -tioj, which are not available in the dc limit. 

Admittedly, the frequency range for resonant enhance- 
ment is difficult to access experimentally. Current mea- 
surements on nanoscopic conductors hardly reach the GHz 
regime |33,34|. Nevertheless, appropriate gating of the 
device could move the HOMO/LUMCGclose to the Fermi 
energy, resulting in resonance enhancement at lower fre- 
quencies. Small gap materials like Graphene nanoribbons 
would offer another route for the experimental realization 
of this effect. 



quencies indeed quadratic and linear scaling, respectively. 

For larger frequencies, resonances appear around 2.5 
eV/h and 4 eW/h with conductances that are two orders of 
magnitude larger than the dc one. This can be qualitatively 
understood by a comparison to the analytical result of Fu 
and Dudley for a single-level model [.24 J . characterized by 
a Breit-Wigner transmission 



TiE) 



7 



{E - EoY + 7^ ' 



(9) 



where £^0 denotes the level energy. Here the admittance is 
given by 



Re{y(o.)} = 



Go 



2uj 



arctan 



f AE + uj 






7 



arctan 



f AE 



(10) 



and 



4 Summary In this study, we investigated the behaviour 
of a contacted 1,4-benzenediol molecule subjected to an 
alternating bias directly and using a Fourier transform of 
Lorentzian and exponential voltage signals. The Lorentzian 
input signal led to very good agreement with reference 
discrete frequency calculations. In the admittance, capaci- 
tive behaviour could be identified and interpreted through 
the transmission of the junction. An analytical single level 
model showed large qualitative similarities to our numeri- 
cal results. The approach we employed for these findings is 
a combination of the highly efficient tight-binding approx- 
imation to adiabatic TDDFT and a device density matrix 
propagation scheme derived within the Keldysh formalism 
by Zheng and co-workers. Following earlier discussions 
on this topic, we can also confirm that a different choice of 
the self-consistent Hamiltonian does not change the equal- 
ity of the TD steady state with its static counterpart at the 
NEGF level. 
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Im{r(w)} = 
^ 7 ^^ / [{^E + L,f + 7^] [{AE uf + 7^] 
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with AE = Ep - Eq. In Fig. |6]b) the Fu-Dudley ad- 
mittance is shown for AE = 2 eV and 7 — 0.25 eV. The 
qualitative similarity to the atomistic results for a real junc- 
tion is clearly seen. Our calculations include a large num- 
ber of molecular states, while Eqs. 10 and 11 hold only 
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